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DESIGN OF RECURSIVE DIGITAL FILTERS HAVING SPECIFIED 
PHASE AND MAGNITUDE CHARACTERISTICS 


By Robert E. King 
Langley Research Center 

and 

Gregory W. Condon 

Langley Directorate, U.S. Army Air Mobility R&D Laboratory 

SUMMARY 

A method for a computer-aided design of a class of optimum filters, having specifi- 
cations in the frequency domain of both magnitude and phase, is described. The method, 
an extension to the work of Steiglitz, uses the Fletcher -Powell algorithm to minimize a 
weighted squared magnitude and phase criterion. Results using the algorithm for the 
design of filters having specified phase as well as specified magnitude and phase compro- 
mise are presented. 


INTRODUCTION 

Recursive filters, wherein the output sequence is both a function of the input as well 
as past output samples, are commonly used in digital signal processing and analysis. 

Such digital filters in many applications offer distinct advantages of precision and versa- 
tility over their continuous or analog counterparts. There exist a number of design pro- 
cedures for implementing digital filters (see ref. 1) each one of which strives to attain 
some analogy between discrete and continuous systems. Transform methods such as the 
matched-z, bilinear -z, and standard-z which lead to specific property invariances are 
available (see ref. 2) to the designer familiar with continuous filter design. 

For frequency-domain synthesis (see refs. 3 and 4), realization is normally by 
means of cascade or parallel combinations of pole and zero pairs in the complex plane. 
The synthesis problem is, in fact, reduced to one of approximation since the filter topol- 
ogy is generally specified. In none of the available design procedures, which can yield 
filters having excellent magnitude -frequency characteristics, however, do the resultant 
filters, in themselves, have particularly useful phase characteristics. Indeed, in striving 
for particular magnitude characteristics by using any of the available design methods, 
there is no control over the filter phase properties. 



In practice, it is often desirable to specify a digital filter in the frequency domain 
by its phase (see ref. 5) or even a compromise between magnitude and phase. The pro- 
cedure in this paper meets these requirements through the use of an iterative computer- 
aided design leading to an optimum set of parameters for a specified filter topology and 
is an extension of the technique described by Steiglitz (see ref. 6) for determining the 
optimum coefficients of a cascade filter having magnitude specifications alone. The 
extension makes possible the design of a new class of digital filters having the prescribed 
phase characteristics. 


SYMBOLS 


A 

4 

< 

4 

e k 

f V 8A 

f k 

f s 

H(z) 

l H k| 

H k 

a | H k|/ 9 P 

K ) 

i,.. .,N 


filter multiplier 

denominator of ith stage of H(z) at £2^ 
magnitude error at £2^ 
phase error at £2^ 
error vector at £2^ 

derivative of error vector at £2^ with respect to zero frequency gain 

frequency at kth specification point, Hz 

sampling frequency, Hz 

unity gain discrete transfer function 

magnitude of H(z) at £2^ 

conjugate of H(z) at £2^ 

gradient vector of magnitude of H(z) at £2^ with respect to parameter 
vector 

imaginary part of quantity 
denotes filter stage 


2 



I 


Jk 

k 

M k 

4 

P 

Pi 

qj(k) 
q^k) 
R( ) 


u i(k) 


V 

v k 

V 

av/aA 

av k /9e k 


w 


k 





w'(k) 


Jacobian at £2^., 


l dp ! 9p_ 


sample point 

specification magnitude at £2^ 
numerator of ith stage of H(z) at £2^ 
parameter vector 

set of filter parameters for the ith stage, a^, b^, c^, and di 

first system state of ith stage at kth sample point 

second system state of ith stage at kth sample point 

real part of quantity 

input to ith stage at kth sample point 

criterion functional, that is, V(A,p) 

criterion functional at £2^., that is, V^(A,p) 

reduced criterion functional, that is, v(A*,p) 

slope of criterion functional with respect to zero frequency gain 

gradient vector of criterion functional at £2^ with respect to error vector 
at £2k 

weighting matrix at £2^- 
magnitude weighting at £2^. 
phase weighting at £2^ 

dummy variable of ith stage at kth sample point 
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Y(z) 

yj(k) 

z 

z k 

0 k 

x 


digital filter discrete transfer function 
output of ith stage at kth sample point 
transform variable 

discrete transform variable at £2^., e 
specification phase at £2jj, radians 
collective phase weight 


0k phase of H(z) at £2^, radians 

80k/ 9 P gradient vector of phase of H(z) at £2^ with respect to parameter vector 
£2^ fractional frequency at kth specification point 


An asterisk on a symbol denotes an optimum value. A circumflex denotes optimi- 
zation with respect to A. A superscript T denotes the transpose. 

DISCUSSION 


The Filter Form 

The fundamental advantages of the N-stage cascade canonical form of recursive 
digital filter whose signal flow graph is shown in figure 1 and which is described by the 
product operator 


N 

Y(z) = A jj 
i=l 


1 + ajZ"l + bjZ ^ 
1 + Cjz"l + djZ - ^ 


( 1 ) 


Y(z) = AH(z) 


are (1) its relative insensitivity to perturbations in the denominator coefficients, an 
important consideration in digital filters, especially of high order and particularly where 
finite register lengths (see ref. 1) are involved; (2) its simplicity of implementation; and 
(3) the simplicity of factoring the filter operator to determine its roots. This form has 
found extensive application in practical filters for signal processing, and a version 
employing serial arithmetic (ref. 7) is commercially available. 
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For completeness, an alternative description of the filter is given in terms of the 
system states and q?j and clearly demonstrates the recursive nature of the filter. 
The set of difference equations describing the filter and required in developing a com- 
puter algorithm is presented. Thus, for the ith stage in figure 1 at the kth sample point 

w i (k) = A i u i (k) - qq^k) - c^q^k) 
q^(k + 1) = w 1 (k) 
q^(k + i) = q^(k) 

Yi(k) = w x (k) + a i q^(k) + b^k) 

where 


u i(k) = yi_i(k) 


is the input to the ith stage and is identical to the output of the (i - 1) stage and 


Ai = 



(i - 1) 


(i^D 


The Synthesis Problem 

The design problem considered in this paper can be stated as follows: When the 
magnitude and phase specifications (M^ and 6^, respectively) at the kth fractional 
Nyquist frequencies $2^ = 2fy i /f s (where f s is the sampling frequency in Hz) are known, 
determine the set of optimum parameters p* of an N-stage cascade filter having the 
form of equation (1) so that the resultant digital filter will have a minimum sum squared 
magnitude and phase error for all specified frequencies. 

By constraining the filter topology, the optimum synthesis problem becomes one of 
parametric optimization with respect to a given criterion of fit. The composite criterion 
which can weight the magnitude and phase requirements independently and as functions of 
frequency is chosen as the inner product 


v (A,p)- ^(e k ,W k e k )= ^ V k 
k k 


( 2 ) 
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where 


e k = 


A ) H k| " M k 


v M 

E k 

^k “ 0 k 




is the error vector and 


W k 



0 


0 

xwf 


is the diagonal weighting matrix. Clearly, V (A,p) is a nonlinear function of the param- 
eter vector p = ^a^,bj,ci,dj,. . ., a N’ b N> c N>dN) T j which involves the 4N filter coeffi- 
cients, and of the filter multiplier A. 


The Minimization Algorithm 

Through formal differentiation of the criterion function (eq. (2)) with respect to the 
multiplier A, the minimization procedure can be slightly simplified to that of finding the 
minimum of a reduced functional V(p) 

Thus 

ay _ v /S 9 V k\ _,Y 
(A { 

and 9V/dA = 0 implies 

| H k| w “( A *l H k| -M k )=0 
k 

or 

^ K|wfM k 

A * = % — r (3) 

2 ki 

k 

An additional necessary condition for existence of an extremum is that the gradient vector 
be zero; thereby, the optimum parameter vector p* is obtained. From equation (2) 


= V (A*,p) involving only 4N parameters. 


Klwf 


e k 
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( 4 ) 


^ = 2Y(j k ,W k e k ) 
a P k x 7 

where the (4N x 2) Jacobian J k is 



8p 


T 


( 5 ) 


Clearly, each element of the gradient vector is the sum of two weighted functions of the 
magnitude and phase error. By writing 


H k| 2 = H k H k 


where H k is the conjugate of H k evaluated at the fractional frequency f2 k , it is 
readily shown (see ref. 6), where p^ is the set of filter parameters for the ith stage, 
that 


8 M _j_ 

“Pi ~| H k| 



For the cascaded filter topology in terms of the elements of p., 




and 
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jirfi k 

where, with z k = e , 

N k = Nl ( z k) = 1 + a i z k 1 + biz k 2 

and 

D k = D X z k) = 1 + c i z k + ^k 2 

By letting 


H k = H k e 


j0 k 


it follows that 


*k ‘ K los e H k) 


whence 


9^k 

dp 


“ '( J- ‘»Se H k 




which takes on a particularly simple form for the cascade topology. For the ith stage 
parameters, in fact, 
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and 


9 ^ k = _! 


9d i \d£ 



The special case of a one -stage (N = 1) filter is illustrated. Here 


H k = A 


1 + az^* + bz^ 
1 + cz' 1 + dz^ 


V = 2 (a* |H k | - M k ) 2 wf + X £ (> k - e k ) 2 w£ 


and 


av 9 

eT = 2 l r k w k -si 

k 


e m w m lil + XE 0 W * ffk N 


"8^- 


/ z -l\ 

qJ^rMt- + XR,fI 

U/ 


? 4 >J 

l k J 


Similarly, 


av = V 
8b Z, 


8V _ V 
8 C Zy 







8V Y 
8d Zy 


/ z -2\ /z" 2 

qMr X UaR^iR- 
ll) k/ W, 


where 


Q“ = 2E“w“|H k | 


and 


= 2E^W^ 

are the weighted errors. It is obvious that the frequency intervals of the input data 
(specifications) need not be uniform and may, in fact, be intentionally unequal to allow for 
nonuniform frequency weighting. 

Complementary Root Reflection and Stability 

In deriving the frequency response of a discrete operator by letting z^ lie on the 
unit circle r, it is possible to take advantage of a unique property of the discrete trans- 
form pertaining to its magnitude when a root lying outside the unit circle is imaged or 
reflected into the unit circle. It is easy to show that the magnitude of a phasor z - zq, 
where zq is a root of the discrete transform lying outside the unit circle, is equal to 


i III 

l 

1 

N 

O 

II 

N 

O 

z 

z 0 


Since zq has been assumed to be outside the unit circle, 1 /zq must be inside, the 
term |zq| correcting for magnitude changes. Thus, if in the optimization procedure a 
pole should stray outside the unit circle and thereby lead to an unstable filter, root reflec- 
tion guarantees stability with no magnitude change. There is no analogous simple identity 
for the phase of a reflected root. Experience with the procedure has shown that provided 
the design requirements can be met by means of a stable filter, that is, that a feasible 
solution exists, an optimum will indeed be found through repeated application of root 
reflection. 


The Computer Algorithm 

A complete listing of the filter design algorithm, which is an adaptation of the pro- 
gram written by Steiglitz, is given in the appendix. The main program is termed STGZ3 
which calls four principal subroutines: (1) FUNCT performs the functional and gradient 
computation for each iteration as well as putting out the final optimum parameters and 
plots, (2) FLPWL is a Fletcher -Powell conjugate gradient routine, (3) INSIDE computes 
root reflection, and (4) ROOTS determines the poles and zeros of the filter. Single- 
precision arithmetic has been employed. 

When minimization of the functional has been attained in the first pass or the mini- 
mization algorithm has iterated 300 times, a test is made to ascertain that all the roots 
are within the unit circle, a necessary requirement for the poles for stability reasons 
and for the zeros to insure minimum phase. If the design should result in an unstable 
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configuration, the roots are reflected about the unit circle and minimization is resumed 
in a second pass. If a minimum does indeed exist and all the roots then lie within the 
unit circle, the program computes and prints out the frequency response and commences 
plotting. 

Minimization is deemed to be achieved when the absolute difference in functionals 
between successive iterations e = |v new - V 0 id| or the norm of the gradient vector 
falls below preassigned limits. Convergence is generally fast for magnitude or phase 
filters but can be very slow for the case of compromise filters. 

When the design specifications cannot be met after LIM iterations (see appendix), 
the program will stop; this situation indicates that the optimum could not be found and the 
resultant characteristic which may be unusable is plotted. Generally, feasible designs 
have been determined in less than 2000 iterations. 

Minimization of the criterion function does not guarantee determination of a global 
minimum but rather determination of a local minimum. Depending upon the parameter 
vector utilized for initialization of the algorithm computation, different minima may be 
achieved. Experience has shown that stage-by-stage optimization, that is, utilization of 
the ith-stage optimum parameter vector as the initial parameter vector for the (i + 1) 
stage of an N-stage filter, yields lower minimum values of the criterion function than 
does single-pass optimization. 


APPLICATIONS 
Linear -Phase Filter 

This example considers a digital filter having application as a phase discriminator 
with a linear phase characteristic and arbitrary magnitude characteristic and is shown 
in figure 2. In this example all magnitude weights were set to zero and all phase weights 
to unity, the multiplier A being arbitrarily made unity since it has no effect on the 
phase characteristic. 

The phase requirements were 0^ = 1 - 212^ (0 = f2jj = 1)> aa d a two-stage filter was 
specified. When an initial parameter vector p = (0, 0, 0, 0.25, 0, 0, 0, 0)T was used, the 
algorithm converged to the optimum, with e = 10 “4, in 52 iterations and a Control Data 
6600 computer time of 14 seconds. The optimum parameter values computed were to 
four places 


A = 1.0 




aj = 0 

b x = -0.9871 

O 

H * 

11 

o 

dj = 0.0395 

a 2 = 0 

b 2 = -0.9871 ' 

o 

11 

CNJ 

o 

d 2 = -0.0127 


I 



It is interesting to note that the phase requirements were met to within 0.00877 radian 
for approximately 95 percent of the frequency range. 

Constant -Phase Filters 

Two cases were considered to obtain filters having constant phases of -tt/ 2 and 
ir/2 radians over a frequency range 0.3 = % = 0.7. As in the previous case, the form 
of the magnitude characteristic was of no concern; hence, zero magnitude weighting was 
specified. With the same initial parameter state used in the previous example, the 
first case (lag network) optimized in 1673 iterations and 42 seconds to yield a hyper- 
bolic magnitude characteristic and phase errors of less than 0.000377 radian throughout 
the specified band. 

The computed parameters for the lag case were 
A = 1.0 


aj = 0.5580 

b t - -0.1857 

cj = -0.4752 

d x = 0.0363 

a 2 = 0.5580 

b 2 = -0.1857 

c 2 = -0.3712 

d 2 = -0.5686 


The positive phase filter (lead network), however, took only 165 iterations and 
17 seconds to yield the desired phase characteristic with errors nowhere exceeding 
O.OOI77 radian in the specified band. 

The optimum filter parameters for this second case were determined to be 
A = 1.0 

a 1 = -0.4768 bi = -0.1548 cj = 0.5022 d x = -0.1082 

a 2 = -0.4768 b 2 = -0.1548 c 2 = 0.4515 d 2 - -0.2008 

It is noted that for both cases, the phase weights outside the specified band were set 
to zero, and thereby allowed for arbitrary phase in these regions. Figures 3(a) and 3(b) 
show the resultant frequency characteristics for the lag and lead cases, respectively, of 
two-stage filters. The combination of the two filters, although they have antagonistic 
magnitude characteristics, suggests the possibility of a phase -splitting digital network. 
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Limited-Band Constant-Gain Linear-Phase Filter 

The third example demonstrates a compromise design of a digital filter having 
constant-magnitude and linear -phase characteristics, over a limited frequency band, 
typical of phase discriminators. Here, except for X = 0, the specifications were stated 
as 


M k = 


(0.3 ^ n k 2 0.7) 

(Elsewhere) 



( 0.3 S f2 k s 0.7) 

(Elsewhere) 


Equal error and frequency weights were employed and the effects of changes in X are 
shown in figure 4 for a two-stage design. Figure 4(a) shows the case of X = 0, that is, 
a magnitude -only filter being specified, and coincidentally yields the linear -phase -filter 
characteristic derived in the first example. (See fig. 2.) Figures 4(b) and 4(c) show the 
magnitude and phase characteristics for the cases of X = 10 and X = 1000, respectively. 
The increasing weight on phase and resultant degradation in the magnitude characteristic 
are shown. The optimum parameters were 

X = 0: 


A = 0.2063 
a^ = 0.0000 
a 2 = 0.0000 


b t = -1.0000 
b 2 = -1.0000 


X = 10: 


A = 0.3658 

aj = -0.9754 bj = 0.7300 

a 2 = 0.8632 b 2 = 0.5632 


cj = 0.0000 

c 2 = 0.0000 


c x = 0.4529 
c 2 = -0.6119 


dj = 0.1539 
d 2 = 0.1539 


dj = 0.7211 
d 2 = 0.7443 
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X = 1000: 


A = 0.4232 

a! = -1.1739 bj = 0.8489 

a 2 = 1.1739 b 2 = 0.8489 


Cj = 0.7596 dj = 0.6691 

c 2 = -0.7596 d 2 = 0.6691 


Low-Pass Zero-Phase Filter 

The fourth example considers a compromise filter, having two and three stages, 
with specifications that are intentionally conflicting. A filter described by 


6> k = 


fl.O 

(0.0 ^ S2 k < 0.5) 

70.5 

(°k-0.5) 

0.0 

V 

(0.5 < S 1.0) 

/° 

(o.o ^ £ 0.5) 

| Unspecified 

(Elsewhere) 


is specified. 

Figures 5 and 6 show the results for the two- and three -stage designs, respectively, 
with figures 5(a) and 6(a) showing the magnitude -only (X = 0) case. The degradation in the 
magnitude characteristics when greater emphasis is placed on the phase specifications 
is evident in figures 5(b) and 6(b) for X = 10 and in figures 5(c) and 6(c) for X = 1000. 
Comparison of figure 6 with figure 5 demonstrates the improvement brought about by 
increasing the number of stages. The optimum parameters for the two-stage filter were 

X =0: 


A = 0.1196 




a! = 1.0240 

b! = 1.0000 

c 1 = -0.1713 

d x = 0.7676 

a 2 = 1.0240 

b 2 = 1.0000 

c 2 = -0.5324 

d 2 = 0.2286 
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X = 10: 


A = 0.4879 



a 1 = 0.2018 

b t = 0.6684 

c j = 0.3560 

d i = 

: 0.4612 


a 2 = 0.6597 

b 2 = 0.4335 

c 2 = 0.0806 

d 2 = 

0.7671 

X = 1000: 

A = 0.5343 






aj = 0.0205 

bj = 0.7169 

cj = -0.0836 

d l = 

0.6255 


a 2 = 0.6286 

b 2 = 0.7905 

c 2 = 0.2123 

d 2 = 

0.6681 

The optimum parameters for the three-stage filter were 
X = 0: 

A = 0.0510 




aj = 0.8537 

bj = 1.0000 

ex = -0.1068 

di = 

1.0000 


a 2 = 0.8537 

b 2 = 1.0000 

c 2 = -0.4046 

d 2 = 

0.5990 


a 3 = 0.8537 

b 3 = 1.0000 

c 3 = -0.6799 

d 3 = 

0.2069 

X = 10: 

A = 0.5109 






a 1 = 1.3302 

b x = 0.5515 

Cj = -0.1731 

dj = 

0.8097 


a 2 = 0.6844 

b 2 = 0.7157 

c 2 = 1.1880 

d 2 = 

0.5850 


a 3 = -0.0373 

b 3 = 0.7012 

c 3 = 0.3825 

d 3 = 

0.5262 
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A = 1000: 


A = 0.4515 
a! = 1.5107 
a 2 = 0.5825 
a 3 = -0.1663 


bj = 0.5286 
b 2 = 0.7490 
b 3 = 0.7485 


cj = -0.1771 
c 2 = 1.3094 
c 3 = 0.2002 


di = 0.8972 
d 2 = 0.4191 
d 3 = 0.6393 


A three-stage design of this example is used to demonstrate the existence of two distinct 
local minima, dependent upon the initial parameter vector. In the first case, a single- 
pass optimization was accomplished with p = (0, 0, 0, 0.25, 0, 0, 0, 0)^ for the initial 
parameter vector and resulted in the optimum filter shown in figure 6(a). In the second 
case, a stage-by-stage optimization was accomplished by utilizing the optimum parameter 
vector from a two-stage design for the initial parameter vector of a three-stage design 
and resulted in the optimum filter shown in figure 7. Comparison of these results demon- 
strates the existence of two distinct local minima, the stage-by-stage minimization 
yielding superior results. 


CONCLUDING REMARKS 

A method has been developed for a computer-aided design of cascade canonical 
digital filters having prescribed magnitude or phase characteristics or a compromise 
between the two. The method, which uses an unconstrained minimization algorithm, 
allows for arbitrary error and frequency weighting. Representative designs of phase and 
compromise filters have demonstrated the utility of the technique. Although convergence 
is generally fast for magnitude phase filters, it may be slow for the case of compromise 
filters. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., February 17, 1972. 
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APPENDIX 


PROGRAM LISTING 


This appendix contains a program listing written for the Control Data 6600 computer 
at the Langley Research Center, Hampton, Virginia, and is an adaptation of that written by 
Kenneth Steiglitz at Princeton University for the design of specified magnitude -only 


filters. 







PRTGR AM STG73 1 INPUT, OUTPUT , TAPE 5= INPUT • TA?F6 = CUT PL T, PUNCH ) 


0013 

ooooo? 


FXTfRNAL FUNCT 


OC 1 4 

00000? 


DIFFUSION Hi 1 84) . X ( X 6) ,G< 16 ) 


CC1 5 

00000^ 


rCPP.?N/RArf/rt< 1C0I .Y< 100 1 ■ K . PHASED f IOC ) . AL AMOA # FR , w TM< ICO 1 . 


0016 



C W TO( 100 1 » K T Y ° 


0C1 7 

000000 


CCP Vf'N/RAWl / I C ALL tKC ALL *L H 



ooooo? 


CALL CALCCMP 


OCX 9 

000004 


CAL L LFR 1Y 


0020 

ooooos 


WP I T c ( 6 ,* 1 1 


0C2 1 

0000 u 

51 

FrfiP A T ( * INPUT DATA* ) 


OC 7 2 

000011 


M = C 


0C2 3 

0000 1 2 

7C 

P = M f 1 


00? 4 

000014 


PFAn(^,21)Vs(M) ,Y( Kl # PH A S E D ( P ) , w T M ( p ) , w rp ( w > 


C02 5 

000031 

21 

F r P P A T ( 5 F 1 0 • 5 ) 


007 6 

0000?1 


WP IT F(6,?2) W ,U(M) ,Y ( H , PHASED! P ) , UTP ( M I ,'*TP ( p > 


0027 

000051 

?2 

FOR P A T ( * I=*.I?»* W= * F 7.4,* Y= * F 7.4,* P H A S >= D= * ,F 7 . 4 , 


002? 



C. * u T M = *,F7.4 , * ^TP = *,F7.4) 


0029 

000051 


I F l »( Ml . LT . 1 . CC IGCTP30 


C 03 0 

000054 


DU 16 J = 1 ,16 


00 31 

000056 

l c 

X ( J ) = 0 .00 


CO 3 ? 

000061 


Xt 4 )= . ?5 


CC 3 7 

03006? 

qc 

R F A C ( 5,60)L ,L I P • F S T , F P S *HPAX,AlA M rA»FR , K t y p 



000106 

60 

F OP P A T(?I5,5FiC.5,nCJ 


CC 35 

0001 06 


IF(FP .1 T. .001 ) F R = 1 . 


OC 3 6 

00011? 


N = 4 *1. 


00 ? 7 

0001 1? 


I F { F H F , 5 ) ^9 G , F 8 « 


OC ? 8 

000117 

58? 

C DN T T MJF 


003° 

0001 17 


W D ITF!6,61IL*LIP,EST ,FPS»FPAX,fK, A L A ^ C A 



0001 41 

61 

F U k V A T ( * L=*.i2»* L I P= * » I 5 , * E S T = * , F 1 0 . 5 , * FPS=*,F10.5, 





C* hPA X=* , F10 . 5 .* FRFCRANCE=*,F10. 5 # * L APS 3 A = * • F 10 , ^ ) 


0042 

0001 41 


ID ALl =0 


CC4 3 

00014? 

Qh 

KCALL =0 



0001 43 


CALL FLPWLI FUNCT , N » X ,F,G,FST ,FPS, c C * IF?, HI 



000155 


CALL prnTS(N.X) 


0045 

0001 57 


CALL I NS I DF 1 N . X, KFLAD ) 


C046 

00016? 


WPITF(6,26)IER,KFLAG,IC£LL,KCALL 



0001 76 

7 6 

FOPPATt* I FR = * » I 5 • * KF L AG = * ,15 » * I C At L ^ , I c » * KCALL=*»I5) 



0001 76 


IF( KCALL. GT.*CC1 GU TO 



ooo?o? 


l F ! (KFLAG.NF . C .0° . If P. \E.O ) . AND. I C ALL . LF. L I P ) GC TC ^ 



00071? 


CALL Rr?TS(N,X) 


CC 5 0 

000715 


rCA(_L =-10 


CC5 1 

000716 


CALL FUNCT t N»X,F,G,HPAX) 


CC C 2 

00022 2 


GU TO c a 


CC " ? 

000223 

9 c .° 

CALL r ALPLKCo t''.* 00 ?! 


CC f L 

000226 


STOP 


cc ct 

000230 


end 


cr r t' 
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SljHR' _ ‘UTiMF r iJMCT ( A *X •F*G*HNAX) 

0057 

GOOD 10 


0 1 m fj" s i c'i r-wrofl ( 200 » ,ph asfx 1200 ) , amagi ?co) . phfi 100 >. *»hac iroi 

0058 

0^0010 


[’ IMFNS FHN XPL 1 2C0. ) *YPLI 20C ) #CX(2 n C ) .CYC ?00 1 

0C5-9 

on DC in 


im'FNS T;v-: Htl S4) . X (I 6) »H16 ) • YHTC 100) . P ( IOC ) 

0060 

0000 i n 


Cf.MPl Ex 7 OP AR » 7 C » 7 7 0 UP ♦ ? 7 f U P 2 

0061 

0000 !0 


r.LMPl FX 7 ( ICO 1 .TUMI ICO f 4 ) .OFNI 10 c .4 > .0. CP AP . ZCURt ZCU«2 .PNFO 

0062 

OO0010 


CCMPf N/RAU/Wt ICC ) .Y( IOC ) .P , FHASFO ( 1 00 ) • At A m DA , F R , Vs TM ( 1 00 ) . 

0063 



C fcTPfiOO ) .ktyp 

0064 

OC 30 10 


rCKMN/OAWl/ICALl .KC ALL .1 IP 


000010 


f A'FC = 0PPLX(1 .00,0.00 ) 

0C6 6 

000012 


P I = - ; • 14l5n?4-7 cqc 7° 

0C6 7 

OOQO 14 


K = N / 4 

0C6 8 

00 10 1- 


IF( If Al L.NF.O 100 Trim 

0069 

00301 S 


or 10 7 1=1 t M 

0070 

00 00 20 

1 32 

Z(I) = CCXPt C^FLX(0.30 .a*( I )*PI ) ) 

0071 

00 f )4G 

101 

Al =0. CO 

0072 

000( 4 l 


A 7 = c • CO 

CC7 3 

ooor 4? 


pn 43 T = 1 * v 

0074 

00004- 


7CIj<= 7(1) 

0075 

00004^ 


7Cl,P2 =70 JP*7CUP 

0076 

OOQn-S '' 


3= O’K X ( 1 .00 ,0.°C> 

0077 

OOQ 0 n 6 


nr ^ J = 1 ,K 

0C7 8 

00)087 


j4= ( j-n * 4 

0079 

OO0("Sl 


TIjWU ,.U = l.o04-X(J4H )*7C UP fX ( J4 h? ) *7CU* ? 

0080 

nOO; 0? 


HE. MI * Jl = 1.0AfX( J4f- )*ZCIJC *X< J4 f4 ) * z C U p 7 

0081 

0031 ?r 

i -- 

o = o * T t j v ( I » j ) /r> p N i i , j ) 

0082 

0031 47 


CP A P = CJ3JGC3) 

0083 

0001 *.? 


YhT ( I ) = O^OPAP 

0084 

0031 40 


A 2 = A 7 h yht ( I) TM 1 ) 

0085 

0001 6 4 


YHT ( I ) - SOP T ( Yh T ( I ) ) 

0086 

0031 6 7 

43 

Ai r A) f YHT 1 I) *Yl II *taTPI I ) 

0087 

000 2 3 1 


I P ( K T Y 15 . *! E • 0 ) Of TO 

0088 

000202 


,\ = A 1 / A ? 

0C8 9 

0^0204 


r,r Tr 64^ 

0090 

00 0 7 0 ^ 

4 fc 4 

A=’ . 

009 1 

0002 0 7 

6 4 7 

CON T T Mu p 

C09 2 

^0070^ 


•3P 4 1 7 1 = 1 » M 

CC9 3 

00C2 11 


ZZCU^ = 7(1) 

0094 

003M4 


77UJP 2=Z7C’JR*7 l DP 

C09 5 

000271 


70=C V,3 LX( 1 . ,3 . ) 

0096 

000274 


PHA S= 0 . 

CC97 

00377=; 


nr 4ii j-i * k 

0098 

OOQ 7 2 <*■> 


J47 =( J-l ) *4 

0099 

000773 


7 07 AL = r\jP O + X ( J 4 7 + 1 >*7ZCUP+X(J^7+? ) *7/ClJP? 

0100 

000 7 <+“’ 


1 Al =70PAu 

0101 

0007 


7 A7=f VP( X ( 0. . - • . ) *70hA p 

Cl 02 

O 03 ' 6 1 


P H A S = P h A S * A T A N 2 ( 7. A 2 . 7 A 1 ) 

0103 

00G*>6 4 


70 = 70*7 )o A p 

0104 

000 7 7 7 


7Q»UR=f N-O + X ( J4Z4-3 l*77CLP + X ( J47+4 ) «/7CUP2 

0105 

OOQ 7 ] ^ 


7A1 =7 3 1 A P 

0106 

r- 

c . 
c 
c 


7 A*5 = r.M^LX( 0. . -1. ) */06AP 

0107 

000377 


PHAS=PHAS-ATAN7(/A2. 7A1 ) 

0108 

000 


7 u = 7 0 / 7 0 H A F 

0109 

0^3147 

4 1 1 

CONTINUE 

0110 

HO0353 


p HA ( I ) = - n H A S / F I 

0111 

00035 ^ 

412 

Phf ( I l=PHA ( ( l-PHASFH I ) 

0112 

000 3 M 


or 4^ J = 1 ♦ I 6 

0113 * 

000361 

- 7 

r, ( J ) = 1 • 33 

0114 

000-3 6 4 


V I = 0 . 

0115 

003 36 c 


V 2 = 0 . 

0116 

000 7 4*-. 


r n 4 ? 1 = \ f v 

0117 

0 30 3 4 


7UjR = 7 ( r ) 

0118 

0003 7 ^ 


7CUP7 = 7C')p*7CUP 

0119 

030777 


YHT I I ) = A * YHT l I ) 

0120 

00040*5 


c { i }=y hT( I )-v ( n 

0121 

000^3 -1 


V 3. = V 1 + E i l )*E( T )*WTM( I) 

012? 

00341 0 


V2=V7 hPHF( I ) * PH F ( I ) *w T F ( I ) 

0123 

00 , )M 4 


p =F (I) *hTM| J ) 

0124 

00341 0 


T = OHF ( T)^WTP( I ) *? . * A 1 . AN* C A 

0125 
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00942? 


nr 4? 7 j = i ,* 

0126 

00 V+ ? 3 


J4 = 1 .1-11 *4 

0127 

000475 


T=2.*R*VhT( I } / TUM ( I , j) 

0128 

000441 


•31 J4 + 1 )=GI J4+: )+ J*ZCL« 

0120 

000451 


r,{ J4 4- ? j = r, c j4 ) to*7f lr? 

0130 

000461 


0 = - *’-‘*YHT ( I » /rCM( l . J ) 

0131 

00047- 


n( ,|4f U = G l J ** 1 1 It )*7f.UR 

0132 

000505 


r.< j4i-4 ) = o( j 4 1 4 ) to*;cup? 

0133 

000*14 

477 

C( MI M'J C 

0134 

00051 ~ 


01 4? J - 1 , K 

0135 

000* 7 0 


.14= ( J-1 ) *4 

0136 

009*?? 


G 1 J4 + 1) =r, ( J4 + 1I *T*AI MAt, 1 7 CLP/ TON ( I , J ) ) 

0137 

ooo* n 7 


(*,( |4t 7 1 = Gt J4 «- ? J A I *AG ( 7 CUP?/ UiM ( t » J H 

0138 

0005 6 a 


ri(j4+i) = ',ij4 + -) + T *A I M0(-; C i;p / i:FMi I ,j)i 

0139 

0^0* 71 


(-7 CLP? 70 5 M I , J 1 ) 

C 1 4 0 

00160* 

-r 7 

f. iMINlJr 

0141 

000613 


F = V 1 * Al_AMDA*V2 

014? 

000616 


I C All = ICALL + 1 

014? 

000620 


KCALL = KCALL H 


000671 


I F ( KC ALL.GT.?CC»FFTURN 


000674 


IF I (IT ALL 710) *10.EO. IC ALL-1 | WR ITF (6,2= ) ICALL , F, < ', ( J) , J = 1 V NI 

0144 

000657 

75 

FLRWATt* CALL NC . * , l 4 , * F=*,F15.8/( 4F!5.R)> 

014* 

000657 


IF( ICALL. GT.20CO ) GO TC 45C 


0036 6'* 


I F ( IC ALL ♦ GT .0 JPETURN 

0147 

0^0665 


GO Tq 449 


000 5 66 

4*C 

I F « ICALL. GT.L IBH 1 GC T f 44C 


000473 


rftur n 



C • » • • 

. . PR in r HUT 

0148 

000673 

444 

wo I TF t 4 ,501 F 

0149 

000701 

50 

FORMAT ( * FINAL FlJNCT If N VALLF =*.F15.R> 

0150 

003701 


WR I TF ( 6 *511 A 

01 c 1 

000707 

51 

FORMATS A=*. F15.fi) 

C 1 r 2 

00 )707 


WP I TF 1 6,5? ) I X ( J ) , J = 1 » N ) 

0153 

000731 

52 

FORMATS FINAL X = */(* *.4F15. fall 

Cl 54 

000731 


WRITF (6.54) (GIJ) . J = 1 ,N ) 

0155 

000751 

54 

FORMAT!* FINAL GRADIENT = */(* «,^F1 C .8M 

0156 

000753 


OH 55 1 = 1. M 

0157 

000760 


WR I TF ( 6.5611 . W ( I) .Y( I 1 . YHT ( I) , E( I ) .PHASED! I ) . PH A 1 I ), PHE! I ) 

0159 

001011 

56 

FORMAT!* I=*,I?,* W=*.FG.5,* Y=*,F«.5,* YHT=* t pp.5 t * E=*,Ffl.5, 

C 15 9 


f. * PHASFO=*,F8.5 Pt-A = *,FF.5,* * F E = * , F « . 5 ) 

0160 

00 10 u 


WR I TF ( 4 . 5 9 ) K 

0161 

001016 

65 

FORMAT!* FINAL TABLE FOR A * . 1 2 . * STAGE 5ILTFP*) 

0162 

001014 


YMIN33=0. 

016 * 

001017 


YHA X3 i = 0. 

0164 

001070 


YMAX1 1=0. 

0165 

001071 


S = F R / 2 00 . 

C 1 6 6 

001073 


nn 60 i=i,?oi 

0167 

001077 


F R F G= S *F L 0 AT ( 1-1 1 

C 16 8 

001033 


ZCIJR= C*XP! C MPL X(0. CO .FRFC*PI ) ) 

0169 

001040 


ZCUR? =7CUR*7CUR 

C 1 7 0 

001045 


0= CMPLXt 1 •OO.C.OCI 

0171 

00105 ) 


PHASF = 0.00 

017i 

001051 


OP 61 J=1 . K 

0173 

001055 


J4= t J-l » *4 

C 1 7 4 

001057 


GhAR = UNF0tX(J4tl)*ZCLRl-X(J4t?)*ZCUR? 

0175 

001 176 


A 1 = QB AP 

0176 

001100 


A? = CMP LX t 0.00 . -1 .00 ) *QP AR 

0177 

0011 10 


PH A SF = PH A S E t ATAN?(A2.A1 ) 

0178 

0011 14 


Q=C*OP A R 

0179 

001177 


0BAR=0NF3tXU4 t- ) *ZCUk+X< J4+4) *ZCLF? 

0180 

0011 44 


A 1 = C w A R 

C 1 8 1 

001146 


A ?= CM PI x( 0.00 ,-l .00) *OBAR 

019? 

001156 


PH A SF = PHASF- ATAN7IA2.A1) 

0183 

00114? 

61 

O-O/OB AP 

0134 

001 1 7 7 


A1 = C* CCNJGIQ) 

0195 

001207 


Al = A* SORT! A 1 ) 

018 6 

00171? 


PH A SF =-PHA S F / P I 

C 1 8 7 

001714 


TMFGA 1 I ) = F R EO 

0198 

001716 


A M A G ( I I = A1 

0189 

001217 


P F A c F X ( I I = PH A S F 

0190 
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0 01??1 


I F ( PHASE— Y M AX 33) 301 *3C1.4C1 

0191 

001776 


4 C 1 YMAX^3 = PhA5F 

0192 

0017 30 


3 C 1 CONTINUE 

0193 

001230 


IF<PHAS C -YPIN33) 402*302*302 

0194 

001 2 33 


4 C 2 YMIN33=PHASF 

0195 

001235 


30? CONTINUE 

0196 

001735 


I F ( A1-YMAX1! ) ?CC*3CC*4C0 

0157 

001240 


4 CO Y M A X 1 1 = A1 

0193 

001242 


3 CO CONTINUF 

0199 

001242 


60 WFITF(6*62)FQEC*PHASF*AI 

0200 

001251 


62 FOR FA T ( * W=*.F15.R** PH A S E / P I = * * F 1 5 . 8 . * YH T = * • F 1 c ' . 6 > 

0201 


C 

RCAl e computations 

0202 

001261 


I YK AX 1 1=YMAX 1 1 f.9999 

0203 

001264 


Y^AX1=HL0AT( I YNAX1 1 ) 

0204 

001266 


TFfHMAX. EG). 0. ) Gfl TO 332 

02C 6 

001767 


YMA XI = HM A X 

0206 

001270 


33‘ 3 CONTINUE 

0207 

001 2 70 


I YKAX33=YMAX3? »*C59 

0203 

001273 


YMAX^=FLDAT { I YFAX-3) 

0209 

001775 


I Y M T N 3 3 = Y M I N ? 3 - . 9 5 9 9 

0219 

001300 


YMI N3=FL1AT ( I YMN33 ) 

0211 

001301 


A X F =1 OH F/I FS/7 ) 

0212 

00 1 3 0 * 


A Y F =1 OHMAGNITUCES 

C213 

0C 1304 


A FR =7 00 . *FR 

0214 

001306 


NFR =700 

0215 


r 

FAGNTTUHE (CHMPUTFO) - PNECUFNCY PLOT 

021 6 

001 30 7 


CAL L I NFOPL T ( 6 *NF P • f FFO A • 1 . A MAG ,l.C. # FR,0.*YNAX!*0.3*-i0*AX^,-10* 

021 7 



C A Y F , Q ) 

0218 


c 

FAGMTUCF (DFSIRFC) - FPFOUFNfV PLOT 

0 219 

001327 


CALI I NFOPLT ( 1 *P . ^s* 1 f Y* 1 . C . . F P , C . . Y*' A X ! , 9 . 5 , -1 C • A XM , - 1 0 , A YM , 1 U 

C2 2 0 

0013^7 


A YM = 1 OH PHASE/FT 

0221 


c 

PHASE-FRFCUENCY PLOT 

0222 

001 151 


CALL I \Fi)PL T ( C ,NP'P . )*'FGA,l*FhASex*l f O.*FP f YKIN3,YFflX3,0.5r-10» 

0223 



C A X K ' . - 1 j . A Y f- , 0 ) 

0224 


c 

PH F { C r 9 T P ^ 0 ) - FREQUENCY p L F T 

0225 

001^7^ 


CALL I NFO PL T ( 1 * •< , V * 1 * F H A ^FO , 1 * r . , F R » Y F I N 3 * Y F AX 3 * 0 . 5 , - 1 0 ♦ 

0226 



C A X F , -19* AY w .11) 

0227 

00141 0 


- f TUt N 

0228 

On lull 


r Mi 

0229 



5lJH ■[ U r I J r f l ^VsL ( F'iNCT . N ♦ X , F , r 

..PST.FPS.L MIT *rFP»H ) 

0230 

00 0-3 1 ^ 

, : I FP 5 TP\i H ( 1 ) * X ( 

n , ct 1 ) 


0231 

oooo i ~ 

C c; M * i N / < A k 1 / ! r A L L 

*KC AL l • L I F 



0000 1 3 

C M L F UNO T( N , X *F . 

g.hnax ) 


0233 

QOOn ?7 

I r ( <C AL L . GT • 3 r o ) 

GO TO 773 



000 n, 7 

IF { IC Al L.GT.L I w ) 

( 0 jr 7 ; 3 



00.30 3 3 

r,^ in 




no in 75 

771 1^0=3 




00 10 37 

fTIJh N 




ooou 

917 1 r R =o 




000 04 ^ 

KHI\T =0 



0235 

00 JO 4 1 

\l 2 = \ + N 



0236 

00 IP 4 ^ 

\ * = \? ► N 



0237 

09004 * 

* * 1 =n 3 ¥ : 



C238 

ono° 45 

: k = n * i 



0239 

000047 

nr 4 j = 1 • n 



0240 

nr-o J ^ ) 

H ( K ) = 1 0 10 



0241 

01015 « 

N J - A - J 



0242 

000^54 

I r ( N.l ) 5.5,7 



0243 

00 10 5-> 

9'*' 3 L = 1 , N J 



0244 

00006 ) 

<l = Kf L 



0245 
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000061 

3 

H(Kl) =0.00 

000066 

4 

< = KL«-i 

00007? 

5 

KT'UNT = <f!UNT + 1 

00 )C 74 


WR I Ti- l n.501 f 

00010: 

5C 1 

P'VVlTl* KHUN T = * . 1 5 ) 

0001 01 


OU0F=F 

OO01 06 


O'J <= | = 1 • N 

0001 07 


K = N + J 

noon n 


H < K ) =• G { J 1 

0001 1 4 


K = K 4 K 

0001 16 


H ( K ) = X ( J ) 

0001 2 1 


K= 1 f N J 3 

0001 


T =c .on 

OOO! ?•» 


DO 3 t = I , M 

00 }1 24 


T--T-MI ) * H ( K ) 

OOO 1 31 


IF( L- J) 0.7,7 

OOO 1 36 

{' 

K = < + 1 

00^1 3 7 


-/I T.7 M 

OCDl 37 

7 

K = K + 1 

000161 


cruiNU" 

0001 46 

9 

H( J )= T 

000150 


0Y=C. 00 

000150 


HNR M= 0 • 00 

000151 


GNR y = 0 .00 

000153 


on 10 J = 1 » N 

0001 54 


HNR M= HNRM+ ABS (H( J) ) 

0001 60 


GNRM = GNRM »- A P S ( G ( J) 1 

0001 63 

10 

0 Y = DY ♦ H ( J ) *GC J ) 

000173 


I F I OY 1 11,51.51 

0001 74 

11 

IF< HNPM/GNRM-FPSI 51 ,51 .12 

000200 

12 

F Y = F 

000201 


AL F A= 2 .00*1 FST-F)/DY 

000204 


AMR HA = 1 .00 

000205 


IF t ALFA! 15.15.13 

0002 C 7 

13 

I F ( ALFA-AMBDA) 14,15 ,15 

000212 

14 

AM B CA = AL F A 

000214 

15 

ALF A= C . 00 

000215 

16 

F X= FY 

000216 


DX=DY 

000 2 2 9 


DC 17 1=1, N 

00022? 

17 

X ( I ) = X II ) +AMRCA*H( I i 

000231 


CALL FUNCTCN , X ,F .G.HNAX > 

000243 


IF( KC ALL.GT.30C) GO TO 724 

000246 


I F C ICALL.GT.L IM GO TO 724 

000251 


GO TO 918 

000251 

724 

I F R =3 

000253 


RFTURN 

000253 

918 

F Y = F 

000254 


OY=C. 00 

000255 


00 18 1 = 1. N 

000257 

18 

DY= CY +G( I )*H( I ) 

000266 


TF(OY) 19,^6.22 

000267 

19 

IFI FY-FX) 20 ,22,2? 

000272 

20 

ambda=ambca»-alfa 

000274 


ALFA= AMBOA 

000275 


FRR0R=1 .610 

000276 


IFI HNR M*AMBDA- ERROR 116.16.21 

000302 

21 

I ER =2 

000304 


RETURN 

000304 

22 

T =0.0 0 

000305 

23 

I F ( AMBOA) 24,36,24 

000306 

24 

7 = 3 .00*1 FX-FY ) / AMBOA +DX + OY 

000314 


ALFA=AMAX1( A B S ( 7 ) t ABS(OX), 

000326 


DALFA=Z/ALFA 


0246 

0247 

0248 

0249 

0250 

0251 

0252 

0253 

0254 

0255 

0256 

0257 

0258 

0259 

0260 
0261 
0262 

0263 

0264 

0265 

0266 

0267 

0268 

0269 

0270 

0271 

0272 

0273 

0274 

0275 

0276 

0277 
027 8 
0279 
0260 
0281 
0282 

0283 

0284 

0285 

0286 
C289 


0292 

0293 

0294 

0297 

0298 
C29 9 

0300 

0301 

0302 

0303 

0304 

0305 

0306 

0307 

ABS(DY) > 0308 

0309 
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000377 


OALFA=LALFA*DALFA-DX/ALF A*CY/ALFA 

0310 

000333 


IF(OALFA) 51,25,25 

0311 

000335 

25 

W=ALFA* SQPTI CALFA) 

0312 

000340 


ALFA= (DY+W-Z » * AM BO A / < D Y f 7 . CO*W-D X ) 

0313 

000351 


On 26 1=1, N 

0314 

000356 

26 

X ( I )= X ( I ) MT-ALFAI*H ( I ) 

0315 

000366 


CALL FUNCT(N.X,F,G,H*AX) 

0320 

000400 


IF(KC ALL.GT.300) GO TH 725 


000403 


TF( ICALL.GT.L IM| GO TO 725 


000406 


GO TO 919 


000406 

725 

I F P =3 


000410 


RFTLP N 


000410 

919 

IF(F-FX) 27,27,28 


00Q4-1 3 

27 

IF(F-FY) 36,36,28 

0323 

0004 16 

28 

OAL FA =0,00 

0324 

000417 


DO 29 T=1,N 

0325 

000421 

29 

DAL FA =DAL F A + G t I)*H( [ ) 

0326 

000430 


I F ( CALFA) 30,33,33 

0329 

000431 

3 C 

I F ( F-rXl 32.31 • 7 3 

0330 

000434 

31 

I F ( DX-C ALFA! 32,36,3? 

0331 

000436 

3? 

F X = F 

0332 

000437 


DX = DA L F A 

0333 

000440 


T = A LF A 

0334 

00044? 


A MB DA = AL F A 

0335 

000443 


GO TO 23 

0336 

000443 

33 

IF(FY-F) 35,34,35 

0337 

00044^ 

34 

I F 1 OY-DALFA ) 35.36,3 5 

0338 

000447 

35 

FY - F 

0339 

000450 


OY = DA L F A 

0340 

000451 


A MR F A = AMBDA-ALFA 

0341 

000454 


GO TO 22 

0342 

000454 

36 

On 37 J = 1 .N 

0343 

000456 


K = N + J 

0344 

000457 


HIK)=Gf Jl-HIK) 

0345 

000463 


K = K fN 

0346 

000464 

37 

HtK ) = X ( 1 ) — H C K ) 

0347 

00047? 


IF(OLDF-F+FPS ) 51,38,38 

0348 

000476 

38 

I FR =0 

0349 

000477 


IF1K0UNT-N) 42,39.39 

0350 

000501 

39 

T = 0 .00 

0351 

00050? 


Z -0 .00 

0352 

0005 O'* 


DC 40 J = 1 , N 

035? 

000504 


K - N f J 

0354 

COO 5 0 C 


rf=H (K ) 

0355 

0C0 C -10 


K = K + \ 

0356 

000*11 


T= T ► APS( H ( K ) 1 

C357 

000514 

40 

7=7trt*H(K) 

0358 

0005? 5 


I F 1 HNJR V-PPS) 4 1,41,42 

0359 

000 K ?6 

41 

TF( T-FPS) 56,56,4? 

0360 

000 c 31 

*T ? 

i r c KruNr-LiMi n 43 , 57 , 5 c 

0361 

00 0 C 34 

43 

ALF A=C .CO 

0362 

000535 


no 47 j = 1 , N 

0363 

0005 *7 


K = J +M 

0364 

000540 


w=c .00 

0365 

00054 ? 


OC 44 L = l ,N 

0366 

000 c 4 7 


Kl =\H. 

0367 

000544 


W = lft+H ( KL ) *HC K 1 

0368 

000 5 5 9 


IF(L-J) 44,45,45 

0369 

0005*4 

44 

K=K f M-L 

0370 

C00557 


GC Tn 46 

0371 

000557 

45 

K=K f 1 

0372 

000561 

46 

CONTI N1LJP 

0373 

000*64 


K=N f J 

0374 

000565 


ALFA=ALFA+W*H(K) 

0375 
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APPENDIX - Continued 


000 c 7 7 
00057ft 
COOftOO 
00060? 
COOft 05 
000605 
OOOftOft 
000607 
000ft?5 
0006 55 
000634 
000ft 3ft 
0006 56 
000640 
000ft 41 
000647 
OOOftft! 
00066 ^ 
000667 
000667 
000671 
0006 71 
000674 
000 A 7ft 
000700 
000700 
000701 
00070? 


000006 

000006 

000006 

000007 

000011 

000014 

000016 

000020 

000022 

000024 

000025 

000027 

000031 

000033 

000035 

000051 

000051 

000054 

000060 

000064 

000066 

000066 

000071 

000071 

000072 

000074 

000076 

000077 


47 

H t J )- W 

0376 


I F ( ALFA ) 48 • 1 .4R 

0377 

4R 

K=N3l 

0378 


00 49 L = 1 • N 

0379 


K L = N2 + L 

0380 


OH 49 J = L • N 

0381 


M J=N?* J 

0382 


H(K)=MK)+HtKL)*HlNJ)/Z~H(L)*HlJ)/ALFA 

0383 

49 

K-K +1 

0384 


GO T 1 ft 

0385 

5 C 

I FP = l 

0386 


P F T UR N 

0387 

51 

on 5 0 J = 1 »N 

0388 


K = N?«- J 

0389 

52 

X ( J ) = H ( K ) 

0390 


r At L FUNC T ( N • X »F»G?HRAX) 
IFI KC ALL.GT.30C) GO TC 726 
I F ( IC ALL.GT.L IM GO TO 7?6 
GO TO ^20 

0393 

7 ? t 

I FR =3 
PFTUP N 


S?C 

IHGMRM-FPSI cc .55.5? 


j 3 

I F C IFR) 56.*4 t ?4 

0396 

6 6 

I F R =- 1 

0397 


Gf’ T^ 1 

0398 

55 

I F P =0 

0399 

56 

RF TUP N 

0400 


FND 

0401 


SUBROUTINE INSIDE (N»X*KFLAG) 

DIMENSION X I 1 6) 

J=-l 
KFL AG=0 
10 J = J+2 

I Ft J. GT • N ) RETURN 
E 5 00 + X ( J ) 

C=X( J + l) 

DISC=B*B-C 

I F C DI SC.LE.O. 00 ) G0T020 

C REAL ROOTS 

D I SC= SQRT(DISC) 

P 1 = B + D I SC 
P 2 = B-DI SC 
CRT = ABS(Rl) 

DR2= ABS(R2) 

IF (ORIOLE. 1.00. AND. DR2.LE. 1.00) GOTO 10 
KFL AG=1 

IF(DR1.GT.1.00)R1=1 .OO/Rl 
IF(DR2.GT.1.00)R2=1 .00/R2 
X(J)=-1.00*<Rl+R2> 

X ( J + l ) =R1*R2 
GCT010 

C COMPLEX ROOTS 

20 IFIC.LE.l .00 ) GO TO 10 
KFLAG=1 
C=1.00/C 
X( J+l )=C 
X< J)=X( J)*C 
GOT 01 0 
END 
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APPENDIX - Concluded 




SUBROUTINE ROOTS(N.X) 

0433 

00000 5 


DIMENSION X ! 1 6 ) 

04 34 

00000 5 


WR lTE( 6 »40 ) 

0435 

000010 

40 

FORMAT ( * R00TS*/6X,*REAL*» 11 X, * I MAG* ,1 1 X, *RE AL*t 11X, *IMAG* » 

0436 

000010 


J^l 

0437 

000011 

10 

J = J + 2 

0438 

0000 1 3 


IF( J.GT.NI RETURN 

0439 

000017 


B=-.500*X( J) 

C440 

000021 


C=X( J+l) 

C441 

000023 


DI$OB*8-C 

0442 

000025 


IF (DISC. LE. 0.00) GO T020 

0443 


c # # Q ^ 

►•REAL RCOTS 

H444 

000027 


DI$C= S ORT (DISC) 

G445 

000030 


Rl-B+DISC 

C446 

00003 2 


R2*B-DISC 

0447 

000035 


WRITE(6»30)R1»R2 

0448 

000044 

30 

FORMAT !* *.F15.8.15X,F15.8) 

0449 

000044 


G0T010 

C450 


c # mm # m 

•COMPLEX ROOTS 


000046 

20 

OISC= SQRT!-1.00*DISC) 

0452 

000052 


DI SCM=-1.00*DISC 

0453 

000055 


WRITE! 6.50) B. DISC. B.DISCM 

0454 

000070 

50 

FORMAT!* *,4F15.8) 

0455 

000070 


G0T010 

0456 

00007 2 


END 

0457 
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(a) Unspecified phase filter. A = 0. 

Figure 4.- Two-stage limited-band constant-gain filters. 
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(b) Linear -phase filter. X = 10. 
Figure 4.- Continued. 
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(c) Linear -phase filter. A = 1000. 
Figure 4.- Concluded. 
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(c) Zero-phase filter. A = 1000. 
Figure 5.- Concluded. 
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Figure 6.- Three-stage low-pass filters. 
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